Local intraspecific aggregation in phytoplankton model communities: spatial scales of occurrence and implications for coexistence.

The coexistence of multiple phytoplankton species despite their reliance on similar resources is often explained with mean-field models assuming mixed populations. In reality, observations of phytoplankton indicate spatial aggregation at all scales, including at the scale of a few individuals. Local spatial aggregation can hinder competitive exclusion since individuals then interact mostly with other individuals of their own species, rather than competitors from different species. To evaluate how microscale spatial aggregation might explain phytoplankton diversity maintenance, an individual-based, multispecies representation of cells in a hydrodynamic environment is required. We formulate a three-dimensional and multispecies individual-based model of phytoplankton population dynamics at the Kolmogorov scale. The model is studied through both simulations and the derivation of spatial moment equations, in connection with point process theory. The spatial moment equations show a good match between theory and simulations. We parameterized the model based on phytoplankters' ecological and physical characteristics, for both large and small phytoplankton. Defining a zone of potential interactions as the overlap between nutrient depletion volumes, we show that local species composition-within the range of possible interactions-depends on the size class of phytoplankton. In small phytoplankton, individuals remain in mostly monospecific clusters. Spatial structure therefore favours intra- over inter-specific interactions for small phytoplankton, contributing to coexistence. Large phytoplankton cell neighbourhoods appear more mixed. Although some small-scale self-organizing spatial structure remains and could influence coexistence mechanisms, other factors may need to be explored to explain diversity maintenance in large phytoplankton.

Flow in temporally and spatially varying porous media: a model for transport of interstitial fluid in the brain.

Flow in a porous medium can be driven by the deformations of the boundaries of the porous domain. Such boundary deformations locally change the volume fraction accessible by the fluid, creating non-uniform porosity and permeability throughout the medium. In this work, we construct a deformation-driven porous medium transport model with spatially and temporally varying porosity and permeability that are dependent on the boundary deformations imposed on the medium. We use this model to study the transport of interstitial fluid along the basement membranes in the arterial walls of the brain. The basement membrane is modeled as a deforming annular porous channel with the compressible pore space filled with an incompressible, Newtonian fluid. The role of a forward propagating peristaltic heart pulse wave and a reverse smooth muscle contraction wave on the flow within the basement membranes is investigated. Our results identify combinations of wave amplitudes that can induce either forward or reverse transport along these transport pathways in the brain. The magnitude and direction of fluid transport predicted by our model can help in understanding the clearance of fluids and solutes along the Intramural Periarterial Drainage route and the pathology of cerebral amyloid angiopathy.

Modelling Plasmid-Mediated Horizontal Gene Transfer in Biofilms.

In this study, we present a mathematical model for plasmid spread in a growing biofilm, formulated as a nonlocal system of partial differential equations in a 1-D free boundary domain. Plasmids are mobile genetic elements able to transfer to different phylotypes, posing a global health problem when they carry antibiotic resistance factors. We model gene transfer regulation influenced by nearby potential receptors to account for recipient-sensing. We also introduce a promotion function to account for trace metal effects on conjugation, based on literature data. The model qualitatively matches experimental results, showing that contaminants like toxic metals and antibiotics promote plasmid persistence by favoring plasmid carriers and stimulating conjugation. Even at higher contaminant concentrations inhibiting conjugation, plasmid spread persists by strongly inhibiting plasmid-free cells. The model also replicates higher plasmid density in biofilm's most active regions.

A Genuinely Hybrid, Multiscale 3D Cancer Invasion and Metastasis Modelling Framework.

We introduce in this paper substantial enhancements to a previously proposed hybrid multiscale cancer invasion modelling framework to better reflect the biological reality and dynamics of cancer. These model updates contribute to a more accurate representation of cancer dynamics, they provide deeper insights and enhance our predictive capabilities. Key updates include the integration of porous medium-like diffusion for the evolution of Epithelial-like Cancer Cells and other essential cellular constituents of the system, more realistic modelling of Epithelial-Mesenchymal Transition and Mesenchymal-Epithelial Transition models with the inclusion of Transforming Growth Factor beta within the tumour microenvironment, and the introduction of Compound Poisson Process in the Stochastic Differential Equations that describe the migration behaviour of the Mesenchymal-like Cancer Cells. Another innovative feature of the model is its extension into a multi-organ metastatic framework. This framework connects various organs through a circulatory network, enabling the study of how cancer cells spread to secondary sites.

An immuno-epidemiological model with waning immunity after infection or vaccination.

In epidemics, waning immunity is common after infection or vaccination of individuals. Immunity levels are highly heterogeneous and dynamic. This work presents an immuno-epidemiological model that captures the fundamental dynamic features of immunity acquisition and wane after infection or vaccination and analyzes mathematically its dynamical properties. The model consists of a system of first order partial differential equations, involving nonlinear integral terms and different transfer velocities. Structurally, the equation may be interpreted as a Fokker-Planck equation for a piecewise deterministic process. However, unlike the usual models, our equation involves nonlocal effects, representing the infectivity of the whole environment. This, together with the presence of different transfer velocities, makes the proved existence of a solution novel and nontrivial. In addition, the asymptotic behavior of the model is analyzed based on the obtained qualitative properties of the solution. An optimal control problem with objective function including the total number of deaths and costs of vaccination is explored. Numerical results describe the dynamic relationship between contact rates and optimal solutions. The approach can contribute to the understanding of the dynamics of immune responses at population level and may guide public health policies.

Analysis of long transients and detection of early warning signals of extinction in a class of predator-prey models exhibiting bistable behavior.

In this paper, we develop a method of analyzing long transient dynamics in a class of predator-prey models with two species of predators competing explicitly for their common prey, where the prey evolves on a faster timescale than the predators. In a parameter regime near a singular zero-Hopf bifurcation of the coexistence equilibrium state, we assume that the system under study exhibits bistability between a periodic attractor that bifurcates from the singular Hopf point and another attractor, which could be a periodic attractor or a point attractor, such that the invariant manifolds of the coexistence equilibrium point play central roles in organizing the dynamics. To find whether a solution that starts in a vicinity of the coexistence equilibrium approaches the periodic attractor or the other attractor, we reduce the equations to a suitable normal form, and examine the basin boundary near the singular Hopf point. A key component of our study includes an analysis of the long transient dynamics, characterized by their rapid oscillations with a slow variation in amplitude, by applying a moving average technique. We obtain a set of necessary and sufficient conditions on the initial values of a solution near the coexistence equilibrium to determine whether it lies in the basin of attraction of the periodic attractor. As a result of our analysis, we devise a method of identifying early warning signals, significantly in advance, of a future crisis that could lead to extinction of one of the predators. The analysis is applied to the predator-prey model considered in Sadhu (Discrete Contin Dyn Syst B 26:5251-5279, 2021) and we find that our theory is in good agreement with the numerical simulations carried out for this model.

A Mathematical Model of TCR-T Cell Therapy for Cervical Cancer.

Engineered T cell receptor (TCR)-expressing T (TCR-T) cells are intended to drive strong anti-tumor responses upon recognition of the specific cancer antigen, resulting in rapid expansion in the number of TCR-T cells and enhanced cytotoxic functions, causing cancer cell death. However, although TCR-T cell therapy against cancers has shown promising results, it remains difficult to predict which patients will benefit from such therapy. We develop a mathematical model to identify mechanisms associated with an insufficient response in a mouse cancer model. We consider a dynamical system that follows the population of cancer cells, effector TCR-T cells, regulatory T cells (Tregs), and "non-cancer-killing" TCR-T cells. We demonstrate that the majority of TCR-T cells within the tumor are "non-cancer-killing" TCR-T cells, such as exhausted cells, which contribute little or no direct cytotoxicity in the tumor microenvironment (TME). We also establish two important factors influencing tumor regression: the reversal of the immunosuppressive TME following depletion of Tregs, and the increased number of effector TCR-T cells with antitumor activity. Using mathematical modeling, we show that certain parameters, such as increasing the cytotoxicity of effector TCR-T cells and modifying the number of TCR-T cells, play important roles in determining outcomes.

Bounds on the Ultrasensitivity of Biochemical Reaction Cascades.

The ultrasensitivity of a dose response function can be quantifiably defined using the generalized Hill coefficient of the function. We examined an upper bound for the Hill coefficient of the composition of two functions, namely the product of their individual Hill coefficients. We proved that this upper bound holds for compositions of Hill functions, and that there are instances of counterexamples that exist for more general sigmoidal functions. Additionally, we tested computationally other types of sigmoidal functions, such as the logistic and inverse trigonometric functions, and we provided computational evidence that in these cases the inequality also holds. We show that in large generality there is a limit to how ultrasensitive the composition of two functions can be, which has applications to understanding signaling cascades in biochemical reactions.

Artificial Neural Network Prediction of COVID-19 Daily Infection Count.

This study addresses COVID-19 testing as a nonlinear sampling problem, aiming to uncover the dependence of the true infection count in the population on COVID-19 testing metrics such as testing volume and positivity rates. Employing an artificial neural network, we explore the relationship among daily confirmed case counts, testing data, population statistics, and the actual daily case count. The trained artificial neural network undergoes testing in in-sample, out-of-sample, and several hypothetical scenarios. A substantial focus of this paper lies in the estimation of the daily true case count, which serves as the output set of our training process. To achieve this, we implement a regularized backcasting technique that utilize death counts and the infection fatality ratio (IFR), as the death statistics and serological surveys (providing the IFR) as more reliable COVID-19 data sources. Addressing the impact of factors such as age distribution, vaccination, and emerging variants on the IFR time series is a pivotal aspect of our analysis. We expect our study to enhance our understanding of the genuine implications of the COVID-19 pandemic, subsequently benefiting mitigation strategies.

Coevolution of Age-Structured Tolerance and Virulence.

Hosts can evolve a variety of defences against parasitism, including resistance (which prevents or reduces the spread of infection) and tolerance (which protects against virulence). Some organisms have evolved different levels of tolerance at different life-stages, which is likely to be the result of coevolution with pathogens, and yet it is currently unclear how coevolution drives patterns of age-specific tolerance. Here, we use a model of tolerance-virulence coevolution to investigate how age structure influences coevolutionary dynamics. Specifically, we explore how coevolution unfolds when tolerance and virulence (disease-induced mortality) are age-specific compared to when these traits are uniform across the host lifespan. We find that coevolutionary cycling is relatively common when host tolerance is age-specific, but cycling does not occur when tolerance is the same across all ages. We also find that age-structured tolerance can lead to selection for higher virulence in shorter-lived than in longer-lived hosts, whereas non-age-structured tolerance always leads virulence to increase with host lifespan. Our findings therefore suggest that age structure can have substantial qualitative impacts on host-pathogen coevolution.

A Mathematical Model for the Impact of 3HP and Social Programme Implementation on the Incidence and Mortality of Tuberculosis: Study in Brazil.

In this paper, we presented a mathematical model for tuberculosis with treatment for latent tuberculosis cases and incorporated social implementations based on the impact they will have on tuberculosis incidence, cure, and recovery. We incorporated two variables containing the accumulated deaths and active cases into the model in order to study the incidence and mortality rate per year with the data reported by the model. Our objective is to study the impact of social program implementations and therapies on latent tuberculosis in particular the use of once-weekly isoniazid-rifapentine for 12 weeks (3HP). The computational experimentation was performed with data from Brazil and for model calibration, we used the Markov Chain Monte Carlo method (MCMC) with a Bayesian approach. We studied the effect of increasing the coverage of social programs, the Bolsa Familia Programme (BFP) and the Family Health Strategy (FHS) and the implementation of the 3HP as a substitution therapy for two rates of diagnosis and treatment of latent at 1% and 5%. Based of the data obtained by the model in the period 2023-2035, the FHS reported better results than BFP in the case of social implementations and 3HP with a higher rate of diagnosis and treatment of latent in the reduction of incidence and mortality rate and in cases and deaths avoided. With the objective of linking the social and biomedical implementations, we constructed two different scenarios with the rate of diagnosis and treatment. We verified with results reported by the model that with the social implementations studied and the 3HP with the highest rate of diagnosis and treatment of latent, the best results were obtained in comparison with the other independent and joint implementations. A reduction of the incidence by 36.54% with respect to the model with the current strategies and coverage was achieved, and a greater number of cases and deaths from tuberculosis was avoided.

Minimal Mechanisms of Microtubule Length Regulation in Living Cells.

The microtubule cytoskeleton is responsible for sustained, long-range intracellular transport of mRNAs, proteins, and organelles in neurons. Neuronal microtubules must be stable enough to ensure reliable transport, but they also undergo dynamic instability, as their plus and minus ends continuously switch between growth and shrinking. This process allows for continuous rebuilding of the cytoskeleton and for flexibility in injury settings. Motivated by in vivo experimental data on microtubule behavior in Drosophila neurons, we propose a mathematical model of dendritic microtubule dynamics, with a focus on understanding microtubule length, velocity, and state-duration distributions. We find that limitations on microtubule growth phases are needed for realistic dynamics, but the type of limiting mechanism leads to qualitatively different responses to plausible experimental perturbations. We therefore propose and investigate two minimally-complex length-limiting factors: limitation due to resource (tubulin) constraints and limitation due to catastrophe of large-length microtubules. We combine simulations of a detailed stochastic model with steady-state analysis of a mean-field ordinary differential equations model to map out qualitatively distinct parameter regimes. This provides a basis for predicting changes in microtubule dynamics, tubulin allocation, and the turnover rate of tubulin within microtubules in different experimental environments.

Integrating Diversity, Equity, and Inclusion into Preclinical, Clinical, and Public Health Mathematical Models.

Mathematical modelling applied to preclinical, clinical, and public health research is critical for our understanding of a multitude of biological principles. Biology is fundamentally heterogeneous, and mathematical modelling must meet the challenge of variability head on to ensure the principles of diversity, equity, and inclusion (DEI) are integrated into quantitative analyses. Here we provide a follow-up perspective on the DEI plenary session held at the 2023 Society for Mathematical Biology Annual Meeting to discuss key issues for the increased integration of DEI in mathematical modelling in biology.

Mathematical Models of Early Hepatitis B Virus Dynamics in Humanized Mice.

Analyzing the impact of the adaptive immune response during acute hepatitis B virus (HBV) infection is essential for understanding disease progression and control. Here we developed mathematical models of HBV infection which either lack terms for adaptive immune responses, or assume adaptive immune responses in the form of cytolytic immune killing, non-cytolytic immune cure, or non-cytolytic-mediated block of viral production. We validated the model that does not include immune responses against temporal serum hepatitis B DNA (sHBV) and temporal serum hepatitis B surface-antigen (HBsAg) experimental data from mice engrafted with human hepatocytes (HEP). Moreover, we validated the immune models against sHBV and HBsAg experimental data from mice engrafted with HEP and human immune system (HEP/HIS). As expected, the model that does not include adaptive immune responses matches the observed high sHBV and HBsAg concentrations in all HEP mice. By contrast, while all immune response models predict reduction in sHBV and HBsAg concentrations in HEP/HIS mice, the Akaike Information Criterion cannot discriminate between non-cytolytic cure (resulting in a class of cells refractory to reinfection) and antiviral block functions (of up to 99 % viral production 1-3 weeks following peak viral load). We can, however, reject cytolytic killing, as it can only match the sHBV and HBsAg data when we predict unrealistic levels of hepatocyte loss.

Mathematical Modelling of Parasite Dynamics: A Stochastic Simulation-Based Approach and Parameter Estimation via Modified Sequential-Type Approximate Bayesian Computation.

The development of mathematical models for studying newly emerging and re-emerging infectious diseases has gained momentum due to global events. The gyrodactylid-fish system, like many host-parasite systems, serves as a valuable resource for ecological, evolutionary, and epidemiological investigations owing to its ease of experimental manipulation and long-term monitoring. Although this system has an existing individual-based model, it falls short in capturing information about species-specific microhabitat preferences and other biological details for different Gyrodactylus strains across diverse fish populations. This current study introduces a new individual-based stochastic simulation model that uses a hybrid τ -leaping algorithm to incorporate this essential data, enhancing our understanding of the complexity of the gyrodactylid-fish system. We compare the infection dynamics of three gyrodactylid strains across three host populations. A modified sequential-type approximate Bayesian computation (ABC) method, based on sequential Monte Carlo and sequential importance sampling, is developed. Additionally, we establish two penalised local-linear regression methods (based on L1 and L2 regularisations) for ABC post-processing analysis to fit our model using existing empirical data. With the support of experimental data and the fitted mathematical model, we address open biological questions for the first time and propose directions for future studies on the gyrodactylid-fish system. The adaptability of the mathematical model extends beyond the gyrodactylid-fish system to other host-parasite systems. Furthermore, the modified ABC methodologies provide efficient calibration for other multi-parameter models characterised by a large set of correlated or independent summary statistics.

Dynamics of Information Flow and Task Allocation of Social Insect Colonies: Impacts of Spatial Interactions and Task Switching.

Models of social interaction dynamics have been powerful tools for understanding the efficiency of information spread and the robustness of task allocation in social insect colonies. How workers spatially distribute within the colony, or spatial heterogeneity degree (SHD), plays a vital role in contact dynamics, influencing information spread and task allocation. We used agent-based models to explore factors affecting spatial heterogeneity and information flow, including the number of task groups, variation in spatial arrangements, and levels of task switching, to study: (1) the impact of multiple task groups on SHD, contact dynamics, and information spread, and (2) the impact of task switching on SHD and contact dynamics. Both models show a strong linear relationship between the dynamics of SHD and contact dynamics, which exists for different initial conditions. The multiple-task-group model without task switching reveals the impacts of the number and spatial arrangements of task locations on information transmission. The task-switching model allows task-switching with a probability through contact between individuals. The model indicates that the task-switching mechanism enables a dynamical state of task-related spatial fidelity at the individual level. This spatial fidelity can assist the colony in redistributing their workforce, with consequent effects on the dynamics of spatial heterogeneity degree. The spatial fidelity of a task group is the proportion of workers who perform that task and have preferential walking styles toward their task location. Our analysis shows that the task switching rate between two tasks is an exponentially decreasing function of the spatial fidelity and contact rate. Higher spatial fidelity leads to more agents aggregating to task location, reducing contact between groups, thus making task switching more difficult. Our results provide important insights into the mechanisms that generate spatial heterogeneity and deepen our understanding of how spatial heterogeneity impacts task allocation, social interaction, and information spread.

Understanding and Quantifying Network Robustness to Stochastic Inputs.

A variety of biomedical systems are modeled by networks of deterministic differential equations with stochastic inputs. In some cases, the network output is remarkably constant despite a randomly fluctuating input. In the context of biochemistry and cell biology, chemical reaction networks and multistage processes with this property are called robust. Similarly, the notion of a forgiving drug in pharmacology is a medication that maintains therapeutic effect despite lapses in patient adherence to the prescribed regimen. What makes a network robust to stochastic noise? This question is challenging due to the many network parameters (size, topology, rate constants) and many types of noisy inputs. In this paper, we propose a summary statistic to describe the robustness of a network of linear differential equations (i.e. a first-order mass-action system). This statistic is the variance of a certain random walk passage time on the network. This statistic can be quickly computed on a modern computer, even for complex networks with thousands of nodes. Furthermore, we use this statistic to prove theorems about how certain network motifs increase robustness. Importantly, our analysis provides intuition for why a network is or is not robust to noise. We illustrate our results on thousands of randomly generated networks with a variety of stochastic inputs.

Advice to a Young Mathematical Biologist.

This paper offers advice to early-mid career researchers in Mathematical Biology from ten past and current Presidents of the Society for Mathematical Biology. The topics covered include deciding if a career in academia is right for you; finding and working with a mentor; building collaborations and working with those from other disciplines; formulating a research question; writing a paper; reviewing papers; networking; writing fellowship or grant proposals; applying for faculty positions; and preparing and giving lectures. While written with mathematical biologists in mind, it is hoped that this paper will be of use to early and mid career researchers across the mathematical, physical and life sciences, as they embark on careers in these disciplines.

Optimal Rotation Age in Fast Growing Plantations: A Dynamical Optimization Problem.

Forest plantations are economically and environmentally relevant, as they play a key role in timber production and carbon capture. It is expected that the future climate change scenario affects forest growth and modify the rotation age for timber production. However, mathematical models on the effect of climate change on the rotation age for timber production remain still limited. We aim to determine the optimal rotation age that maximizes the net economic benefit of timber volume in a negative scenario from the climatic point of view. For this purpose, a bioeconomic optimal control problem was formulated from a system of Ordinary Differential Equations (ODEs) governed by the state variables live biomass volume, intrinsic growth rate, and area affected by fire. Then, four control variables were associated to the system, representing forest management activities, which are felling, thinning, reforestation, and fire prevention. The existence of optimal control solutions was demonstrated, and the solutions of the optimal control problem were also characterized using Pontryagin's Maximum Principle. The solutions of the model were approximated numerically by the Forward-Backward Sweep method. To validate the model, two scenarios were considered: a realistic scenario that represents current forestry activities for the exotic species Pinus radiata D. Don, and a pessimistic scenario, which considers environmental conditions conducive to a higher occurrence of forest fires. The optimal solution that maximizes the net benefit of timber volume consists of a strategy that considers all four control variables simultaneously. For felling and thinning, regardless of the scenario considered, the optimal strategy is to spend on both activities depending on the amount of biomass in the field. Similarly, for reforestation, the optimal strategy is to spend as the forest is harvested. In the case of fire prevention, in the realistic scenario, the optimal strategy consists of reducing the expenses in fire prevention because the incidence of fires is lower, whereas in the pessimistic scenario, the opposite is true. It is concluded that the optimal rotation age that maximizes the net economic benefit of timber volume in P. radiata plantations is 24 and 19 years for the realistic and pessimistic scenarios, respectively. This corroborates that the presence of fires influences the determination of the optimal rotation age, and as a consequence, the net economic benefit.

Untangling the Molecular Interactions Underlying Intracellular Phase Separation Using Combined Global Sensitivity Analyses.

Liquid-liquid phase separation is an intracellular mechanism by which molecules, usually proteins and RNAs, interact and then rapidly demix from the surrounding matrix to form membrane-less compartments necessary for cellular function. Occurring in both the cytoplasm and the nucleus, properties of the resulting droplets depend on a variety of characteristics specific to the molecules involved, such as valency, density, and diffusion within the crowded environment. Capturing these complexities in a biologically relevant model is difficult. To understand the nuanced dynamics between proteins and RNAs as they interact and form droplets, as well as the impact of these interactions on the resulting droplet properties, we turn to sensitivity analysis. In this work, we examine a previously published mathematical model of two RNA species competing for the same protein-binding partner. We use the combined analyses of Morris Method and Sobol' sensitivity analysis to understand the impact of nine molecular parameters, subjected to three different initial conditions, on two observable LLPS outputs: the time of phase separation and the composition of the droplet field. Morris Method is a screening method capable of highlighting the most important parameters impacting a given output, while the variance-based Sobol' analysis can quantify both the importance of a given parameter, as well as the other model parameters it interacts with, to produce the observed phenomena. Combining these two techniques allows Morris Method to identify the most important dynamics and circumvent the large computational expense associated with Sobol', which then provides more nuanced information about parameter relationships. Together, the results of these combined methodologies highlight the complicated protein-RNA relationships underlying both the time of phase separation and the composition of the droplet field. Sobol' sensitivity analysis reveals that observed spatial and temporal dynamics are due, at least in part, to high-level interactions between multiple (3+) parameters. Ultimately, this work discourages using a single measurement to extrapolate the value of any single rate or parameter value, while simultaneously establishing a framework in which to analyze and assess the impact of these small-scale molecular interactions on large-scale droplet properties.